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Thermoelectric  self-cooling  systems  hold  good  prospects  for  the  future,  since  they  improve  the  cooling  of 
any  heat-generating  device  without  electricity  consumption.  The  potential  number  of  applications  seems 
to  be  enormous,  hence  the  necessity  of  a  specific  model  to  simulate  this  type  of  thermoelectric 
application. 

This  paper  presents  a  computational  model  for  thermoelectric  self-cooling  applications,  capable  of 
simulating  both  the  steady  and  the  transient  state  of  the  whole  system.  Supported  by  fluid-dynamics 
software  and  based  on  the  implicit  finite  differences,  the  model  solves  the  system  of  equations 
composed  of  the  Fourier’s  law  and  the  thermoelectric  effects  Seebeck,  Peltier,  Joule  and  Thomson, 
including  temperature-dependant  properties.  Furthermore,  the  model  architecture  allows  the  inclusion 
of  new  analytical  expressions  and/or  procedures  in  order  to  simulate  any  component  with  higher  ac¬ 
curacy  or  include  more  complex  ones. 

Statistical  studies  indicate  ±12%  of  maximum  deviation  between  experimental  and  simulated  values 
of  the  main  outputs.  Furthermore,  the  model  simulates  the  performance  of  the  system  under  abruptly 
changing  conditions. 

In  conclusion,  the  computational  model  turns  out  to  be  a  powerful  tool  that  will  play  a  key  role  in  the 
design  and  development  of  thermoelectric  self-cooling  applications. 

©  2013  Elsevier  Ltd.  All  rights  reserved. 


1.  Introduction 

Thermoelectric  generation  has  been  included  among  the  non¬ 
fossil  fuel  technologies  that  can  play  a  significant  role  in  the 
global  fight  against  greenhouse  effect  and  climate  change,  mainly 
because  it  allows  the  efficient  recovery  of  waste  heat  for  electric 
power  generation  [1,2], 

A  thermoelectric  generator  is  a  solid-state  heat  engine  that  ex¬ 
changes  heat  with  a  cold  and  a  hot  reservoir,  producing  a  certain 
amount  of  work  by  Seebeck  effect.  It  is  a  special  heat  engine  that 
requires  neither  real  fluids  nor  moving  parts,  thus  being  more 
robust,  compact  and  noiseless  than  common  electric  power  gen¬ 
erators  [1,3],  This  fact  explains  why  thermoelectric  generators  are 
working  as  electric  power  suppliers  for  well-known  aerospace 
vehicles,  such  as  Voyager  I  &  II  and  rover  MSL-curiosity,  and  are 
being  installed  in  remote  terrestrial  locations,  such  as  inaccessible 
stretches  of  oil  pipelines,  where  regular  maintenance  is  highly 
difficult  or  even  impossible  [1—3]. 
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However,  the  low  efficiency  of  the  thermoelectric  conversion  is 
acting  as  a  deterrent  to  the  widespread  use  of  this  technology  in  the 
civilian  market  Great  efforts  are  being  made  to  improve  this  effi¬ 
ciency  by  increasing  the  so-called  figure  of  merit  of  thermoelectric 
materials.  However,  researchers  insist  that  further  investigation 
must  be  conducted,  not  only  on  advanced  materials  for  thermo¬ 
electric  conversion,  but  also  on  new  applications  that  make  a  good 
use  of  those  materials  currently  available  on  the  civilian  market  [  1,4], 

In  line  with  this  last  statement,  a  novel  thermoelectric  appli¬ 
cation,  called  thermoelectric  self-cooling  (TSC),  was  presented  in  a 
previous  paper  [5],  This  system  uses  thermoelectric  technology  to 
improve  the  cooling  of  a  heat-generating  device  and  keep  its 
temperature  under  control,  without  electricity  consumption.  It 
works  as  follows:  Given  a  device  that  generates  a  certain  amount  of 
heat  and  is  attached  to  a  passive  cooling  system,  several  thermo¬ 
electric  modules  are  installed  between  them.  These  modules  har¬ 
vest  part  of  the  heat  generated  by  the  device  and  transform  it  into 
electricity,  which  is  supplied  to  the  passive  cooling  system.  As  a 
result,  the  cooling  system  becomes  active,  improving  the  cooling 
power  without  electricity  consumption. 

In  the  cited  paper  [5],  a  TSC  prototype  was  constructed  and 
studied.  Fig.  1  shows  a  sketch  and  a  picture  of  this  prototype.  It 
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ThomsonThomson  effect 

comprises  a  heat-generating  device,  a  finned  dissipator  acting  as 
passive  cooling  system,  four  thermoelectric  modules  that  supply 
power  to  an  axial  fan  installed  over  the  dissipator,  and  a  cold 
extender  that  separates  the  device  from  the  dissipator  to  prevent 
thermal  bridges  between  them.  Experimental  results  point  out 
that  the  thermal  resistance  between  the  heat  source  and  the 
environment  is  reduced  by  30%  compared  to  that  obtained  with 
the  passive  cooling  system,  without  electricity  consumption, 
which  proves  the  success  of  the  application  and  its  good  prospects 
for  the  future. 

The  present  paper  delves  into  this  research  line  with  the  main 
objective  of  developing  a  computational  model  for  simulation  of 
TSC  systems.  The  number  of  potential  applications  seems  to  be  so 
large  that  a  simulation  tool  would  be  of  great  importance  to  assess 
the  potential  impact  of  this  technology  in  the  real  world. 

It  must  be  noted  that,  although  the  final  objective  of  a  TSC 
system  is  the  proper  cooling  of  a  device,  it  comprises  the  common 
components  of  a  thermoelectric  generator,  namely  heat  source, 
heat  sink,  heat  exchangers  and  thermoelectric  modules  acting  as 
electricity  generators.  As  stated  in  several  papers  [6-8],  featuring 
detailed  and  precise  reviews  on  the  state-of-the-art  of  computa¬ 
tional  models  for  thermoelectric  generators,  the  simulation  of  a 
complete  generator  is  a  difficult  task,  given  the  complexity  of  the 
heat  transfer  equations  and  the  influence  of  the  temperature- 
dependant  thermoelectric  effects.  Furthermore,  this  difficulty  in¬ 
creases  when  the  heat  exchangers  must  be  simulated  in  detail,  as 
occurs  in  TSC  applications,  given  their  significant  influence  on  the 
final  performance  of  the  system  [9],  This  is  the  reason  why  most  of 


the  authors  studying  complete  thermoelectric  generators  have 
used  powerful  fluid-dynamics  software  for  the  heat  exchangers  but 
simplified  models  for  the  thermoelectric  modules  —  which  have 
served  their  purpose,  though  -  [7,10-13],  Up  to  now,  only  few  of 
them  have  simulated  the  whole  system  using  fluid-dynamics 
software,  with  notable  success  indeed.  This  last  approach  is  due 
to  predominate  in  the  future,  when  more  powerful  personal  com¬ 
puters  decrease  the  enormous  time  that  the  simulation  requires 
nowadays  [14], 

The  design  and  optimization  of  a  TSC  system  requires  not  only 
detailed  representations  of  the  heat  exchangers  (in  our  case,  the 
cooling  system)  in  both  steady  and  transient  states,  but  also  of  the 
thermoelectric  modules,  including  all  the  thermoelectric  effects 
with  temperature-dependant  properties.  All  our  expertise  in 
simulation  of  thermoelectric  applications  [9,15-19]  indicates  that 
the  combination  of  fluid-dynamics  software  for  the  exchangers  and 
finite-differences  models  for  the  thermoelectric  modules  meets  all 
these  requirements. 

Therefore,  the  present  paper  sets  out  to  develop  a  finite- 
differences  model  for  simulation  of  TSC  systems  in  both  steady 
and  transient  states,  including  all  the  thermoelectric  effects  with 
temperature-dependant  properties,  and  supported  by  fluid- 
dynamics  software  for  the  cooling  system.  Section  2  presents  the 
model,  describing  the  internal  architecture  and  the  representation 
of  every  component  of  a  TSC  system.  Then,  Section  3  shows  the 
verification  and  validation  process,  based  on  statistical  procedures, 
that  provides  the  deviations  of  the  main  output  parameters.  To 
conduct  this  process,  the  values  of  these  outputs  provided  by  the 
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computational  model  are  compared  with  experimental  results 
available  in  the  previous  TSC  paper  [5]. 

2.  Computational  model 

The  finite-differences  method  emerges  as  a  powerful  tool  for 
heat  transfer  simulation,  joining  simplicity,  accuracy  and  high  rate 
of  convergence  [20],  Based  on  this  method,  the  computational 
model  solves  the  non-linear  set  of  equations  composed  of  the 
Fourier’s  law  and  the  thermoelectric  effects  Seebeck,  Peltier,  Joule 
and  Thomson,  all  of  them  temperature-dependant. 

Before  that,  the  TSC  system  must  be  transformed  into  a  mesh  of 
representative  nodes,  each  one  connected  to  the  surrounding 
nodes  by  thermal  resistances  and  endowed  with  a  thermal  capacity 
that  sets  the  variation  rate  of  its  temperature  [20,21],  Fig.  2  shows 
the  mesh  of  27  nodes  that  represents  the  TSC  system,  as  Section  2.3 
explains  in  detail,  whose  performance  is  determined  by  the  Fou¬ 
rier’s  law  and  the  thermoelectric  effects. 

2.1.  Fourier's  law 


The  model  is  based  on  the  implicit  finite-differences  method, 
which  does  not  require  any  restriction  either  on  the  time  step  or  the 
mesh  size  to  converge  [20,21  ].  The  three-dimensional  Fourier’s  law, 
provided  by  Eq.  (1),  is  replaced  with  the  corresponding  expression 
in  implicit  finite-differences.  Eq.  (2)  shows  this  expression  applied 
to  a  node  “a”  connected  to  n  additional  nodes.  Then,  Eqs.  (3)  and  (4) 
provide  respectively  the  thermal  capacity  of  node  “a”  and  the 
thermal  resistances  between  this  node  and  the  n  additional  ones, 
where  U  stands  for  the  corresponding  convective,  conductive  or 
contact  heat  transfer  coefficient.  Once  Eqs.  (2)— (4)  are  applied  to  all 
the  nodes  of  the  system,  a  linear  set  of  equations  arises,  which  can 
be  easily  solved  by  any  numerical  method,  such  as  the  Gaussian 
elimination  or  the  Gauss-Seidel  methodology. 
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It  must  be  noted  that  Eq.  (1)  is  the  special  case  of  the  three- 
dimensional  Fourier’s  law  that  considers  constant  thermal  conduc¬ 
tivity  [21  ].  Derived  from  it,  Eq.  (2)  considers  also  constant  thermal 
conductivity  but  only  inside  each  node  at  a  specific  time  instant,  since 
this  equation  is  applied  locally.  Therefore,  the  thermal  conductivity  is 
allowed  to  vary  between  nodes  and  time  instants,  and  can  be  defined 
as  a  function  of  the  temperature,  as  Section  2.3  explains. 

2.2.  Thermoelectric  effects 

A  thermoelectric  module  comprises  N  thermoelectric  pairs, 
composed  of  p-type  and  n-type  doped  semiconductor  legs,  con¬ 
nected  electrically  in  series  by  metallic  shunts  and  thermally  in 
parallel.  At  the  hot  ends  of  these  pairs,  heat  is  absorbed  by  Peltier 
effect,  whereas  heat  is  generated  at  the  cold  ends.  Eqs.  (5)  and  (6) 
represent  these  effects.  Likewise,  Eq.  (7)  provides  the  heat  gener¬ 
ated  by  Joule  and  Thomson  effects  along  the  2 N  semiconductor  legs, 
being  the  electrical  resistance  of  the  thermoelectric  module  pro¬ 
vided  by  Eq.  (8),  and  the  Thomson  coefficients  derived  from  Seebeck 
coefficients  using  Eq.  (9).  The  electrical  resistance  of  the  module 
comprises  the  electrical  resistances  of  the  legs  plus  the  surface 
electrical  resistances  of  the  soldered  joints  between  legs  and  shunts. 


Q.Peltier,hot  —  NITj,ot  (flip  Ctn)  hot  (5) 

Qpeltier.cold  =  NITco|d(Op  -  ffn)coid  (6) 

Qjoule, Thomson  =  (0  ^0  ~  ^(zp  ~  Tn)  f  (^hot—  ^cold)  (2) 

R0  =  N(ppLp/Ap  +  pnL„/A„  +  2pp/Ap  +  2  psn/An)  (8) 

t  =  T  da/AT  (9) 


Finally,  Eq.  (10)  provides  the  electromotive  force  produced  by  a 
thermoelectric  module,  while  Eq.  (11)  shows  the  electric  power 
generated  when  the  module  is  connected  to  a  specific  load 
resistance. 

E  =  N(((ap  -  an)hotThot  -  (ap  -  an)coldTcold)  -  (tp  -  rn) 

x  (Thot  —  Icold))  (10) 


(4)  P  -  £2Rl/(Rl  +  R0)2 


Ra,i  =  VUajAa  1  = 


(11) 
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Fig.  2.  General  representation  and  electrical  analogy  of  a  TSC  system. 


As  occurred  with  the  thermal  conductivity  in  Section  2.1,  the 
electric  resistivities  and  the  Seebeck  and  Thomson  coefficients  vary 
between  nodes  and  time  instants,  and  can  be  defined  as  functions 
of  the  temperature.  As  a  consequence,  all  the  heat  fluxes  and 
electric  variables  become  temperature-dependant,  and  their  values 
change  along  the  simulation  time. 

2.3.  Simulation  of  TSC  systems 

As  can  be  seen  in  Fig.  1,  a  TSC  system  is  composed  of  a  device 
endowed  with  an  internal  generation  of  heat,  several  thermoelec¬ 
tric  modules,  a  cold  extender,  a  passive  cooling  system  and  a  fan. 
Fig.  2  shows  the  electrical  analogy  of  the  system,  displaying  also  the 
thermal  resistances,  thermal  capacities  and  heat  generation  rates 
relative  to  every  node.  This  electrical  analogy  is  the  gist  of  the 
computational  model  and  is  described  in  detail  along  the  following 
sections. 

2.3.1.  Device  and  heat  source 

In  the  abstract,  any  heat-generating  device  could  include  a  TSC 
system.  However,  given  the  temperature  limitations  of  thermo¬ 
electric  modules  available  in  the  civilian  market,  applications  are 
constrained  to  devices  working  at  low  temperatures.  For  instance, 
electrical  power  converters  and  transformers  are  devices  endowed 
with  internal  heat  generation  that  usually  require  cooling  systems 
in  order  to  keep  their  working  temperature  below  150  °C,  being 
therefore  potential  applications  of  TSC  technology. 

As  can  be  seen  in  Fig.  2,  the  first  two  nodes  on  the  left  represent 
respectively  the  heat  source  and  the  device,  being  Qhs  the  heat  flow 


rate  introduced  by  the  former.  This  heat  is  transferred  from  the  heat 
source  to  the  external  surface  of  the  device  through  the  conductive 
thermal  resistance  Rk, dev-  Then,  part  of  it  is  transferred  to  the 
thermoelectric  modules  through  RCont,hot,  which  stands  for  the 
contact  thermal  resistance  between  the  surface  and  the  thermo¬ 
electric  modules.  The  rest  is  absorbed  by  the  environment,  being 
Rc, dev  the  convective  thermal  resistance  between  the  surface  of  the 
device  and  the  environment.  All  these  thermal  resistances  can  be 
calculated  by  Eq.  (4),  being  the  corresponding  areas  and  heat 
transfer  coefficients  input  parameters  of  the  computational  model. 

Finally,  the  thermal  capacities  of  the  heat  source  (Chs)  and  the 
device  (Cdev)  are  provided  by  Eq.  (3),  being  the  masses  and  specific 
heats  the  last  inputs  that  characterize  the  heat  source  and  the 
device. 

2.3.2.  Thermoelectric  modules 

The  thermoelectric  modules  harvest  part  of  the  heat  generated 
by  the  device  and  transform  it  into  electricity  by  Seebeck  effect. 
Each  module  is  composed  of  several  thermoelectric  pairs  con¬ 
nected  electrically  in  series  by  metallic  shunts  and  thermally  in 
parallel,  being  all  of  them  sandwiched  between  two  plates  that 
provide  electrical  insulation  and  mechanical  robustness.  Any 
thermoelectric  generator  includes  several  of  these  modules,  usually 
connected  also  electrically  in  series  and  thermally  in  parallel. 

2.32.1.  Thermoelectric  pairs.  The  thermoelectric  pairs  are  the 
most-difficult-to-simulate  components,  since  several  effects  occur 
at  the  same  time:  the  four  thermoelectric  effects  Seebeck,  Peltier, 
Joule  and  Thomson,  as  well  as  heat  conduction  based  on  Fourier’s 
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law.  This  is  the  reason  why  the  computational  model  represents 
each  thermoelectric  pair  with  twenty  nodes,  ten  for  each  leg,  as 
Fig.  3  displays.  Each  node  represents  a  ninth  of  a  leg,  except  for  the 
ones  at  either  end,  each  one  of  which  represents  an  eighteenth  of  a 
leg  plus  half  of  a  shunt. 

Since  a  thermoelectric  generator  includes  equal  modules  con¬ 
nected  thermally  in  parallel,  all  the  n-doped  legs  present  similar 
temperature  distributions  and  all  the  p-doped  legs  present  also 
similar  temperature  distributions.  Therefore,  in  the  electrical 
analogy,  all  the  n-doped  legs  can  be  joined  in  one  single  branch 
containing  10  nodes  and,  similarly,  all  the  p-doped  legs  can  be 
represented  by  10  nodes  in  another  branch,  as  Fig.  2  displays.  Then, 
the  corresponding  thermal  resistances  and  thermal  capacities  are 
calculated  by  the  following  equations,  derived  from  Eqs.  (3)  and  (4): 


Rn(i,i+ 1) 

-  L"/9  i  =  J ... 

MNkn(i)A„ 

,9 

(12) 

^p(i,i+l) 

-  W9  t  „  ,  ... 

MN  kp:j)Ap 

,9 

(13) 

Cn(i)  = 

MNA„Lndnc„/9  i  =  2,  - 

•,9 

(14) 

Cp(i)  = 

MNApLpdpCp/9  i  =  2,  • 

■■,9 

(15) 

Cn(i )  = 

MN(A„L„dncn/ 18  +AshL 

sh^shcsh/2)  t  =  1,10 

(16) 

Cp(i)  = 

MN  (ApLpdpCp/1 8  +AshL, 

shdshcsh/2)  1=1,10 

(17) 

Lengths,  cross  areas,  conductivities,  densities  and  specific  heats 
of  legs  and  shunts,  as  well  as  the  number  of  modules  and  pairs  per 
module,  are  inputs  of  the  computational  model.  Moreover,  both 
thermal  conductivities  can  be  defined  as  sixth-order  polynomials  of 
the  temperature. 

Regarding  the  heat  flow  rates,  those  produced  by  Joule  and 
Thomson  effects  are  distributed  along  the  twenty  nodes.  In  addi¬ 
tion,  the  nodes  at  either  end  also  include  the  heat  absorbed  or 
generated  by  Peltier  effect  and  the  heat  generated  by  Joule  effect 
at  the  soldered  joints  between  legs  and  shunts.  Equations  from 
(18)— (23)  provide  the  heat  flow  rates  at  each  node. 


Qnl  =  MAf[  ^  2an')Tn'l  +  ff  +  PnlI2LnH8 


Qpl  =  MN -^)^f  +  ^/2+Ppl/2^g 


fn(i-l)  ~  Jn(i+1)\ 


Qn(i)  =  MN^pn(i)f2L^9 

it 

Q„io  =  MN  [(“pio  -“mo)7nio/  +  J/2  +  Pnlof2 W18 


t  —  2,  ,  9 

(21) 


Q,,o  -  -  «m£mli.il2  +  l,p 


(23) 


Seebeck  coefficients,  electrical  resistivities  and  surface  electrical 
resistivities  of  n-doped  and  p-doped  legs  are  inputs  of  the 
computational  model,  being  the  first  two  ones  defined  as  sixth- 
order  polynomials  of  the  temperature. 

Finally,  the  electrical  resistance  of  the  N  modules  connected 
electrically  in  series  and  the  electromotive  force  generated  by  them 
are  provided  respectively  by  Eqs.  (24)  and  (25),  while  Eqs.  (26)— 
(28)  provide  the  voltage,  electric  current  and  electric  power 
generated  by  the  modules  when  connected  to  a  load  resistance. 

Ro  ■mn[I^(x+t+Sp«‘ i)  - 2X 

+Emo)  +2£]  <«> 


Tp(i-1)-Tp(i+1) 

“  2-s  (JpCO - 1 - an& -  2 - ; )  I 


-  E _ 

Rl  +^o 

E 

Rl  +  Ro 


(25) 

(26) 
(27) 


Fig.  3.  Discretization  of  a  thermoelectric  pair. 


(Rl  +  RoY 


(28) 
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2.3.22.  Electrical  insulation  plates.  The  insulation  plates  at  the  hot 
end  of  the  thermoelectric  modules  are  joined  and  represented  by  a 
single  node,  and  the  same  occurs  with  the  insulation  plates  at  the 
cold  end,  as  Fig.  2  points  out.  Then,  Eqs.  from  (29) — (32)  present  the 
corresponding  thermal  resistances  and  thermal  capacities,  being 
the  length,  cross  area,  conductivity,  density  and  specific  heat  of 
these  insulator  plates  the  last  inputs  of  the  computational  model 
that  define  the  thermoelectric  modules. 


Qns,hot  =  ^ins,hotAns, hotshot  Wot  (31) 


Qns,cold  —  M^ins.coldl-ins.cold^ins.coldCins.cold  (32) 

2.3.3.  Cold  extender 

A  cold  extender  must  be  installed  to  space  the  device  from  the 
dissipator  in  order  to  prevent  thermal  bridges  between  them  that 
would  reduce  both  the  electricity  production  and  the  system  effi¬ 
ciency.  A  simple  metallic  prism  does  this  job  in  the  prototype  dis¬ 
played  in  Fig.  1. 

This  component  is  represented  by  a  single  node  in  the  electrical 
analogy  and  is  connected  to  the  surrounding  nodes  by  RCont,coid. 
which  stands  for  the  contact  thermal  resistance  between  the  mod¬ 
ules  and  the  cold  extender;  Rk,ex,  representing  the  heat  conduction 
through  this  component;  and  Rc  ex,  standing  for  the  heat  convection 
between  the  lateral  surface  of  the  extender  and  the  environment.  All 
these  thermal  resistances  can  be  calculated  by  Eq.  (4). 

Finally,  Eq.  (3)  provides  the  thermal  capacity  of  the  cold 
extender  (Cex).  The  corresponding  dimensions,  heat  transfer  co¬ 
efficients,  density  and  specific  heat  are  inputs  of  the  computational 
model. 

2.3.4.  Cooling  system 

A  plane  plate,  a  finned  dissipator,  heat  pipe,  etc.  could  be  used  as 
passive  cooling  system.  Then,  a  fan  installed  over  it  would  provide 
forced  convection,  so  that  the  cooling  system  would  switch  form 
passive  to  active,  improving  the  cooling  power  without  electricity 
consumption,  since  the  fan  would  be  supplied  directly  by  the 
thermoelectric  modules.  In  the  prototype  presented  in  Fig.  1,  a 
finned  dissipator  acts  as  passive  cooling  system. 

The  dissipator  and  the  fan  are  represented  by  a  single  node  in 
the  electrical  analogy.  The  thermal  resistance  of  this  component 
(i?cS)  includes  the  thermal  contact  between  the  dissipator  and  the 
cold  extender,  the  heat  conduction  through  the  dissipator  and  the 
heat  convection  between  the  dissipator  and  the  environment.  This 
last  term  is  not  obvious,  since  it  depends  on  the  air  flow  forced 
through  the  fins,  which  in  turn  depends  on  the  electric  power  that 
the  modules  supply  to  the  fan.  Therefore,  each  TSC  application  will 
require  a  particular  study  to  determine  the  relation  between  the 
thermal  resistance  of  this  component  and  the  electric  power  sup¬ 
plied  by  the  modules  to  the  fan.  This  relation  is  also  an  input  of  the 
computational  model,  and  therefore  it  must  be  obtained  previously. 
Appendix  A  shows  a  methodology  to  obtain  the  cited  relation  that 
includes  the  use  of  computational-fluid-dynamics  software  and 
experimental  procedures. 

Finally,  the  thermal  capacity  of  the  cooling  system  can  be 
derived  from  Eq.  (3),  being  inputs  of  the  model  the  corresponding 
dimensions,  density  and  specific  heat. 


2.3.5.  Model  architecture 

Fig.  4  presents  the  flow  chart  of  the  computational  model. 
The  input  parameters  are:  dimensions,  densities  and  specific 
heats  of  all  the  components  of  the  system;  heat  flow  rate  pro¬ 
duced  by  the  heat  source;  environment  temperature;  number  of 
thermoelectric  modules  and  number  of  thermoelectric  pairs  per 
module;  thermal  conductivity,  Seebeck  coefficient,  electrical  re¬ 
sistivity  and  surface  electrical  resistivity  of  n-doped  and  p-doped 
legs  as  functions  of  the  temperature;  heat  transfer  coefficients 
(contact,  convection  and  conduction);  electric  load  resistance  of 
the  fan  and  thermal  resistance  of  the  cooling  system  as  a 
function  of  the  electric  power  provided  by  the  modules  to  the 
fan.  Furthermore,  since  the  model  simulates  also  the  transient 
state  of  the  TSC  system,  the  initial  temperatures  of  all  the  nodes 
are  likewise  required. 

The  computational  model  calculates  the  temperature  of  every 
node  for  instant  of  time  t,  updates  the  temperature-dependant 
parameters  and  determines  the  corresponding  thermal  re¬ 
sistances,  thermal  capacities  and  heat  flow  rates.  All  these  pa¬ 
rameters  are  considered  constant  during  a  time  step  At,  and  are 
finally  used  to  calculate  the  new  temperatures  for  instant  of  time 
t  +  At.  The  system  is  considered  to  have  reached  the  steady  state 
when  the  temperature  difference  between  two  consecutive  time 
instants  is  lower  than  a  tolerance  defined  by  the  user.  Every  step 
does  represent  a  real  state  of  the  system,  which  means  that  the 
model  simulates  also  the  transient  state.  It  should  be  noted  that 
the  model  is  deterministic,  since  no  randomness  is  included  in  the 
inputs. 

The  model  outputs  are:  voltage,  electric  current  and  electric 
power  generated  by  the  modules  (and  consumed  by  the  fan);  ef¬ 
ficiency  of  the  modules  and  the  system;  temperatures  of  all  the 
nodes  and  heat  flow  rates  (heat  transferred  between  nodes  and 
generated/absorbed  by  Joule,  Thomson  and  Peltier  effects).  All 
these  outputs  are  presented  as  functions  of  time. 

Finally,  it  is  convenient  to  note  that  the  implicit  finite- 
difference  method  guarantees  the  convergence  of  the  solution 
no  matter  what  time  step  or  mesh  size  is  selected  [20,21], 
Furthermore,  the  electrical  analogy  can  be  easily  extended  to 
obtain  more  refined  meshes,  and  new  analytical  expressions  and 
procedures  can  be  added,  in  order  to  simulate  any  component 
with  higher  accuracy  or  include  more  complex  ones,  as  can  be 
deduced  from  previous  papers  describing  the  use  of  implicit 
finite-difference  methods  to  simulate  thermoelectric  systems 
[9,15-19], 

3.  Verification  and  validation  of  the  computational  model 

The  verification  and  validation  process  (V&V)  assesses  the  ade¬ 
quacy  and  reliability  of  the  computational  model,  that  is,  it  checks 
whether  or  not  the  model  provides  coherent  values  of  the  main 
outputs  and  calculate  the  deviations  from  experimental  values  of 
these  outputs  [22,23], 

All  the  experimental  data  used  to  conduct  the  V&V  process  is 
available  in  the  cited  previous  TSC  paper  [5],  where  one  can  find  the 
design  and  performance  of  the  TSC  prototype  presented  in  Fig.  1.  In 
Section  3.1,  this  prototype  is  described  in  order  to  provide  the 
values  of  the  inputs  that  the  computational  model  requires  to 
simulate  this  TSC  application. 

3.1.  TSC  prototype 

3.1.1.  Thermoelectric  modules 

The  prototype  includes  four  thermoelectric  modules  Kryotherm 
TGM-287-1.0-1.5,  connected  electrically  in  series  and  thermally  in 
parallel.  Each  module  comprises  287  pairs  of  n-doped  and  p-doped 
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Fig.  4. 


bismuth-telluride  legs,  connected  electrically  in  series  by  cupper 
shunts,  and  sandwiched  between  two  similar  aluminium-oxide 
layers.  The  parameters  that  the  computational  model  requires  to 
simulate  these  modules  can  be  seen  in  Table  1  [3,24], 

3.3.2.  Device  and  heat  source 

It  consists  of  an  aluminium  plate,  with  dimensions 
220  x  160  x  22  mm3,  that  encloses  two  thin-film  heating  resistors 


Table  1 

Characteristics  of  the  thermoelectric  module  Kryotherm  TGM-287-1.0-1.5. 


M 

4 

dp  =  d„  6892  kg/m3 

N 

287 

Cins,hot=  850  J/kg  K 

/ti„S.hot  = 

40  x  40  mm 

Cp  ’=  c„  544  J/kg  K 

AP  =  A„ 

lxl  mm 

dtp  =  -an  10  6(-0.002025T2 

+  1. 4234487—44.95361 1 )  V/K 

0.8  mm 

kp  =  k„  0.000029T2— 0.01 9593T  + 

4.809677  W/mK 

Lp  =  L„ 

15  mm 

ft,  =  ft,  10~6(0.043542T— 2.754139)  flm 

dins  hot  = 
dins  cold 

3900  kg/m3 

ps  0.11  flm2 

acting  as  heat  source.  Each  heating  resistor  occupies 
80  x  80  x  0.5  mm3,  providing  up  to  150  W  of  heat  flow  rate  with 
maximum  working  temperature  of  130  °C.  The  density  and  specific 
heat  of  the  aluminium  are  2710  kg/m3  and  890  J/kg  K  respectively, 
whereas  the  main  component  of  the  heating  resistors  is  a  silicon 
layer  with  density  1190  kg/m3  and  specific  heat  1000  J/kg  K.  All 
these  parameters  define  the  thermal  capacities  of  the  device  and 
the  heat  source,  which  turn  out  to  be  1867.8  and  7.6  J/K 
respectively. 

Regarding  the  thermal  resistances,  the  conductive  thermal 
resistance  Rk, dev  was  calculated  experimentally  [5]  and  turned  out 
to  be  0.132  K/W.  Likewise,  the  convective  thermal  resistance  Rk, dev 
is  given  by  Eq.  (33),  where  the  expression  for  the  convective  heat 
transfer  coefficient  was  also  calculated  experimentally  [5],  During 
the  simulation,  the  model  updates  this  parameter  at  every  time 
step.  Finally,  since  all  the  contact  surfaces  between  the  compo¬ 
nents  of  the  system  are  covered  with  conductive  paste  in  order  to 
enhance  the  heat  conduction,  the  contact  heat  transfer  coefficient 
can  be  fixed  to  15,625  W/m2K  [25],  Therefore,  Eq.  (4)  indicates 
that  the  contact  thermal  resistance  between  the  surface  of 
the  device  and  the  four  thermoelectric  modules  RCont,hot  yields 
0.01  K/W. 
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Rcdev  ~  Ac~UcAev 

Acdev  =  Asurf  -  4  x  Ans.hot  =  0.08712  -  4  x  0.0016  =  0.08072  (m* 2) 
Uc4ev  =  6.86  +  0.0369 (lsurf  -  Tenv)  -  0.000056  (Ysurf  -  Tenv)2  (W/m2K) 


(33) 


Finally,  the  heat  flow  rate  generated  by  the  source  (Qhs)  is  used 
as  input  variable  in  the  V&V,  as  can  be  seen  in  Section  3.2. 

3.1.3.  Cold  extender 

It  consists  of  a  50-mm-long  aluminium  prism,  with  base  area 
equal  to  each  module’s  base  area  multiplied  by  the  number  of 
modules  (=0.0064  m2).  The  aluminium  exhibits  thermal  conduc¬ 
tivity  of  202.4  W/mK  and  density  and  specific  heat  similar  to  those 
of  the  device,  indicated  in  Section  3.1.2.  These  values  lead  to  a 
thermal  capacity  (Cex)  of  771.8  J/K  and  a  conductive  thermal 
resistance  {Rk.ex)  of  0.039  K/W,  as  Eqs.  (3)  and  (4)  indicate.  Finally, 
the  convective  thermal  resistance  between  the  lateral  surface  of  the 
extender  and  the  environment  (Rc,ex)  is  calculated  by  Eq.  (33),  being 
now  the  convective  area  equal  to  0.008  m2,  which  corresponds  to 
the  lateral  surface  of  the  cold  extender. 

3.1.4.  Cooling  system  (dissipator  and  fan) 

The  dissipator  is  made  out  of  aluminium,  similar  to  that  used  for 
the  device  and  the  cold  extender,  and  comprises  a  square  base 
plate,  with  side  length  155  mm  and  height  12  mm,  and  23  fins  with 
dimensions  155  x  23  x  1.5  mm3.  A  wind  tunnel  and  a  fan  SUNON 
KDE1208PTS1  cover  the  dissipator.  The  fan  consumes  2.6  W  of 
electrical  power  at  maximum  voltage  (12  V),  and  requires  at  least 
0.4  W  to  start  rotating.  If  the  fan  is  working,  it  exhibits  60  Q  of  load 
resistance;  if  not,  the  load  resistance  falls  to  17  Q. 

The  volume  of  the  dissipator  turns  out  to  be  3  x  10-4  m3,  which 
multiplied  by  the  corresponding  density  and  specific  heat  leads  to  a 
thermal  capacity  (Ccs)  of  723.6  J/K. 

Finally,  as  explained  in  Section  2.3.4,  the  thermal  resistance  of 
the  dissipator  (Rcs)  depends  directly  on  the  electric  power  supplied 
by  the  thermoelectric  modules  to  the  fan.  For  the  present  appli¬ 
cation,  we  devised  a  procedure  that  includes  fluid-dynamics 
simulation  and  experimental  tests  to  obtain  the  relation  between 
RCs  and  the  voltage  supplied  by  the  modules  to  the  fan  ( V ). 
Appendix  A  shows  an  in-depth  explanation  of  this  procedure, 
which  provides  the  relation  indicated  by  Eq.  (34). 


between  the  heat  source  and  the  environment  (ThS— Tenv),  and  the 
temperature  difference  between  ends  of  the  modules  (Tjn- 
s, hot- fins, cold),  all  of  them  measured  once  the  system  reaches  the 
steady  state.  The  experimental  study  conducted  with  the  TSC  pro¬ 
totype,  and  presented  in  the  cited  paper  [5],  provides  experimental 
values  of  these  three  outputs  for  eight  different  values  of  the  heat 
flow  rate  introduced  by  the  heat  source  (Qhs ).  specifically  60, 80, 100, 
120, 130, 160, 190  and  220  W.  These  eight  runs  are  replicated  twice, 
so  that  the  complete  experiment  comprises  twenty-four  runs. 
During  all  of  them,  the  environment  temperature  was  fixed  to  0  °C. 

Now,  the  computational  model  provides  simulated  values  of 
these  three  outputs  for  the  cited  eight  values  of  heat  flow  rate.  As 
mentioned  before,  the  model  is  deterministic,  so  no  variation  in  the 
results  occurs  for  replicated  runs. 

The  V&V  includes  statistical  methods  based  on  confidence  in¬ 
tervals  for  each  one  of  the  three  outputs  [23],  For  example, 
regarding  the  electric  power  supplied  by  the  modules  to  the  fan,  for 
each  heat  flow  rate,  there  are  three  experimental  values  and  one 
simulated.  Then,  if  the  simulated  value  falls  inside  the  confidence 
interval  for  the  mean  formed  with  the  three  experimental  values, 
no  significant  difference  is  considered  between  the  simulated  and 
the  experimental  values.  If  this  occurs  for  the  eight  heat  flow  rates, 
the  model  is  considered  to  predict  accurately  the  electric  power 
supplied  by  the  modules  to  the  fan.  Finally,  if  this  procedure  is 
successfully  applied  also  to  the  other  two  outputs,  the  V&V  of  the 
model  is  complete.  In  terms  of  V&V  criteria,  this  procedure  remains 
as  follows: 

1-  For  each  heat  flow  rate  introduced  by  the  heat  source,  the  99% 
confidence  interval  for  the  mean,  formed  by  the  three  steady- 
state  experimental  values  of  the  electric  power  supplied  by 
the  thermoelectric  modules  to  the  fan,  must  include  the  simu¬ 
lated  value. 

2-  For  each  heat  flow  rate  introduced  by  the  heat  source,  the  99% 
confidence  interval  for  the  mean,  formed  by  the  three  steady- 
state  experimental  values  of  the  temperature  difference  be- 


If  P  >  0.4  W  (Fan  rotating) 

Rcs  =  l/(  _  0.0038V4  +  0.1412V3  -  1.9399V2  +  11.805V  -  22.023)  (K/W) 
If  P  <  =  0.4  W  (Fan  not  rotating) 

Rcs  =  1.2(K/W) 


3.2.  Methodology  and  V&V  criteria 

The  outputs  used  to  conduct  the  V&V  are  the  electric  power 
supplied  by  the  modules  to  the  fan  (P),  the  temperature  difference 


tween  the  heat  source  and  the  environment,  must  include  the 

simulated  value. 

3-  For  each  heat  flow  rate  introduced  by  the  heat  source,  the 

99%  confidence  interval  for  the  mean,  formed  by  the  three 
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steady-state  experimental  values  of  the  temperature  difference 
between  ends  of  the  modules,  must  include  the  simulated 
value. 

Finally,  the  V&V  concludes  with  the  statistical  analysis  of  the 
deviations  between  experimental  and  simulated  values  of  the  three 
outputs,  which  are  calculated  by  Eq.  (35).  This  analysis,  as  well  as 
the  calculations  of  the  confidence  intervals,  is  conducted  with  the 
statistical  software  Statgraphics  Plus  5.1. 

dev  =  100exprSim  (35) 


3.3.  Results  and  analysis 

The  eight  simulations  are  carried  out  with  60  s  of  time  step,  a 
tolerance  of  0.01  K  and  the  initial  temperatures  of  all  the  nodes 
equal  to  the  environment  temperature.  Table  2  shows  the  steady- 
state  experimental  values  [5],  the  simulated  ones  and  the  de¬ 
viations  between  them  for  the  three  considered  outputs,  as  well  as 
the  corresponding  99%  confidence  interval  for  the  mean  (Cl). 

It  can  be  seen  that  all  the  simulated  values  are  included  in  the 
corresponding  confidence  intervals,  indicating  that  the  V&V  of  the 
computational  model  is  achieved.  Therefore,  the  model  is  proven  to 
provide  coherent  values  of  the  outputs,  being  also  accurate  enough 
to  meet  the  three  established  V&V  criteria.  Then,  the  statistical 
analysis  of  the  deviations  provides  the  following  results: 

In  the  first  place,  Statgraphics  indicates  that  the  series  of  de¬ 
viations  between  experimental  and  simulated  values  of  the  tem¬ 
perature  difference  between  the  heat  source  and  the  environment 
in  steady  state  follows  a  normal  distribution  with  mean  /r  =  0.65 


and  standard  deviation  a  —  1.67.  Then,  basic  statistics  [26]  points 
out  that  95%  of  these  deviations  fall  into  the  interval  (-2.62, 3.92), 
as  Eq.  (36)  indicates. 


IC95%  =  (/i  -  1 .96ff,  n  +  1 .96(7)  (36) 

The  temperature  difference  between  the  heat  source  and  the 
environment  is  the  most  important  output  parameter  for  TSC  ap¬ 
plications,  since  the  aim  of  these  systems  is  to  keep  this  tempera¬ 
ture  difference  under  control.  Therefore,  the  computational  model 
stands  out  as  a  powerful  tool  for  the  design  and  optimization  of  TSC 
systems,  since  it  predicts  this  output  with  deviations  included  in  a 
narrow  interval,  specifically  6.54%. 

A  similar  analysis  for  the  electric  power  points  out  that  the  se¬ 
ries  of  deviations  approximates  a  normal  distribution  with 
mean  -1.38  and  standard  deviation  3.07.  Therefore,  95%  of  these 
deviations  fall  in  the  interval  (-7.40,  4.64).  Likewise,  the  series  of 
deviations  of  the  temperature  difference  between  ends  of  the 
modules  follows  a  normal  distribution  with  mean  —0.64  and 
standard  deviation  2.96,  which  means  that  95%  of  these  deviations 
fall  in  the  interval  (-6.44,  5.16).  The  computational  model  predicts 
these  two  outputs  in  steady  state  with  similar  deviations:  12.04% 
for  the  electric  power  and  11.60%  for  the  temperature  difference, 
which  can  be  accepted  given  the  less  significance  of  these  outputs. 
Furthermore,  this  similarity  was  expected,  since  both  outputs  are 
closely  related,  as  Eqs.  (10)  and  (11)  point  out. 

As  mentioned  throughout  Section  2,  the  computational  model  is 
also  capable  of  representing  transient  states  of  the  system.  To  assess 
the  adequacy  and  reliability  of  the  model  in  predicting  the  transient 
performance  of  the  main  outputs,  three  different  tests  were 
proposed: 


Table  2 

Experimental  and  simulated  steady-state  values  of  the  outputs  for  the  V&V  of  the  computational  model. 

Qhs  (W)  Ths-Taw  (K) _  P(W) _  Tins, hot-Wold  (K) _ 

Exp  Sim  Dev  (%)  Exp  Sim  Dev  (%)  Exp  Sim  Dev  (%) 


60  50.4 

50.8  48.9 

49.8 

Cl  (47.4,  53.2) 

80  65.1 

63.4  63.7 

63.5 

Cl  (58.8,  69.5) 

100  78.7 

77.6  77.8 


3.07 

3.89 

1.84 

2.20 

-0.47 


1.16 

-0.26 


120 


130 


190 


Q  (73.1,  82.5) 

91.1 

92.9  91.5 

92.2 

Cl  (86.7,  97.3) 

76.7 

75.5  74.6 

75.4 

CI(71.7,  80.0) 

89.2 

87.8  89.6 

87.9 

Cl  (83.8,  92.8) 

101.9 

100.1  101.9 


-0.44 

1.53 

0.77 

2.82 

1.21 

1.07 


Cl  (95.3, 106.3) 

220  113.8  2.61 

112.6  110.9  1.52 

112.9  1.80 

0(109.3,  116.9) 


0.120  -0.83 

0.122  0.121  0.83 

0.122  0.83 

Cl  (0.115,  0.128) 

0.191  -4.50 

0.186  0.200  -7.00 

0.186  -7.00 

Cl  (0.171,  0.204) 

0.272  -5.56 

0.269  0.288  -6.60 

0.277  -3.82 

Cl  (0.249,  0.296) 

0.389  0.78 

0.374  0.386  -3.11 

0.392  1.55 

Cl  (0.330,  0.440) 

0.702  0.705  -0.43 

0.717  1.70 

Cl  (0.603,  0.798) 

1.132  -0.26 

1.134  1.135  -0.09 

1.134  -0.09 

Cl  (1.127,  1.140) 

1.647  -1.38 

1.640  1.670  -1.80 

1.653  -1.02 

Cl  (1.610, 1.683) 

2.239  2.38 

2.235  2.187  2.19 

2.258  3.25 

a  (2.173,  2.314) 


7.5  7.6 
7.4 

Cl  (7.1,  7.8) 

9.2 

9.7  9.9 

9.6 

Cl  (7.6, 11.3) 

11.9 

12.0  12.1 

11.9 

Cl  (11.6, 12.3) 

14.1 

13.6  14.2 

14.3 

Cl  (11.9, 16.1) 

24.1 

25.0  25.2 

24.4 

0(21.9,27.1) 

30.9 

31.0  29.8 

31.5 

Cl  (29.3,  33.0) 

36.8 

36.9  36.3 

35.9 

Cl  (33.4,  39.7) 

42.7 

42.4  41.7 

42.0 

Cl  (40.4,  44.4) 


2.40 


0.72 


T(K) 

100 
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P(W) 


starting  process  of  the  fan  with  160  W  of  heat  flow  rate  introduced  by  the  heat  source. 

The  first  one  is  presented  in  Fig.  5  and  represents  the  starting 
process  of  the  fan,  showing  experimental  and  simulated  values  of 
Ths- Ten*  Tins, he  >t— Tins, cold  and  P  for  Q,hs  equal  to  160  W.  It  can  be  seen 
that  the  electric  power  rises  as  the  temperature  difference  between 
the  heat  source  and  the  environment  increases.  However,  it  is  not 
until  the  electric  power  reaches  0.4  W  that  the  fan  begins  to  rotate. 
Then,  the  load  resistance  of  the  fan  increases  abruptly  from  17  to 
60  Q  and  so  does  the  electric  power  and  the  temperature  difference 
between  the  ends  of  the  modules.  From  that  moment  on,  forced 
convection  is  obtained  over  the  dissipator,  reducing  its  thermal 
resistance,  which  explains  the  significant  reduction  in  the  slope  of 
the  ThS— Tenv  curve.  Finally,  the  three  outputs  tend  to  stabilization. 

The  second  test  is  presented  in  Fig.  6  and  represents  the  oppo¬ 
site  situation  to  that  of  Fig.  5.  Now,  the  system  is  steady  with  the  fan 
working  and  the  heat  flow  rate  set  at  125  W.  Then,  the  fan  is 
deliberately  blocked.  The  first  consequence  is  that  the  thermal 
resistance  of  the  dissipator  increases  as  natural  convection  takes 
place,  which  explains  the  increase  in  Ths-Tenv-  Likewise,  the  load 
resistance  of  the  fan  falls  to  17  Q  and  the  electric  power  plummets. 
Finally,  when  the  fan  is  released,  the  electric  power  does  not  sur¬ 
pass  the  0.4  W  required  to  operate  the  fan,  though  the  heat  flow 
rate  remains  125  W.  The  three  outputs  stabilize  at  different  values 
from  those  before  the  fan  was  blocked. 


Fig.  6.  Simulated  (broken  curves)  and  experimental  (solid  curves)  values  for  the 
stopping  process  of  the  fan  with  125  W  of  heat  flow  rate  introduced  by  the  heat  source. 


Fig.  7.  Simulated  (broken  curves)  and  experimental  (solid  curves)  values  for  steady- 
state  conditions  with  220, 160, 190  and  again  220  W  of  heat  flow  rate,  and  the  cor¬ 
responding  transitions. 


Finally,  the  third  test  is  presented  in  Fig.  7,  which  shows  a  series 
of  experiments  with  different  heat  flow  rates.  When  the  system  is 
working  steadily,  the  heat  flow  rate  is  changed  abruptly  to  a  new 
value.  Then,  the  system  stabilizes  and  the  process  is  repeated.  This 
figure  shows  four  steady-state  conditions  (heat  flow  rates  of  220, 
130, 160  and  again  220  W)  and  the  corresponding  three  transient 
states. 

All  these  tests  prove  that  the  computational  model  simulates  the 
transient  state  of  the  TSC  system  with  enough  accuracy.  The  simu¬ 
lated  values  of  the  three  outputs  agree  quite  well  with  the  experi¬ 
mental  data,  even  in  abruptly  changing  conditions  such  as  the 
starting  or  stopping  process  of  the  fan.  In  this  sense,  it  is  clear  that 
more  refined  meshes  and/or  smaller  time  steps  would  lead  to  more 
accurate  solutions,  increasing  also  the  simulation  time,  which  would 
not  be  actually  a  problem,  since  every  simulation  is  completed  in 
few  seconds.  However,  the  model  provides  accurate-enough  values 
of  the  main  outputs  in  both  steady  and  transient  state,  which  in¬ 
dicates  that  is  unnecessary  to  conduct  such  refinements  for  the 
studied  configuration  of  TSC  system.  Therefore,  the  present  version 
of  the  computational  model  serves  as  a  design  and  optimization 
tool,  powerful  enough  for  this  type  of  thermoelectric  application. 

4.  Conclusions 

Thermoelectric  self-cooling  systems  hold  good  prospects  for  the 
future  and  are  expected  to  replace  common  forced-convection 
systems  in  the  cooling  of  heat-generating  devices,  since  they 
improve  significantly  the  cooling  power  without  electricity  con¬ 
sumption.  The  potential  number  of  applications  seems  to  be 
enormous,  hence  the  necessity  of  a  computational  model  for  this 
type  of  thermoelectric  application. 

This  paper  presents  a  computational  model  for  TSC  systems, 
capable  of  simulating  both  steady  and  transient  states.  Based  on  the 
implicit  finite  differences,  the  model  solves  the  system  of  equations 
composed  of  the  Fourier’s  law  and  the  Seebeck,  Peltier,  Joule  and 
Thomson  thermoelectric  effects.  Seebeck  coefficients,  electrical 
resistivities  and  thermal  conductivities  of  the  n-doped  and  p-doped 
thermoelectric  legs  are  introduced  as  functions  of  the  temperature. 
Likewise,  we  devised  a  procedure  to  define  the  thermal  resistance 
of  the  dissipator  as  a  function  of  the  voltage  supplied  by  the 
modules  to  the  fan,  which  can  be  used  for  other  TSC  applications. 

The  computational  model  turns  out  to  be  a  powerful  design  tool, 
providing  coherent  values  of  the  most  important  outputs.  Thus,  95% 
of  the  deviations  between  experimental  and  simulated  values  of 


1124 


A  Martinez  et  al.  /  Energy  55  (2013)  1114-1126 


the  temperature  difference  between  the  heat  source  and  the 
environment  fall  into  the  interval  (—2.62%,  3.92%).  Likewise,  the 
interval  (-7.40%,  4.64%)  includes  95%  of  the  deviations  between 
experimental  and  simulated  values  of  the  electric  power  supplied 
by  the  thermoelectric  modules  to  the  fan,  whereas  the  interval 
(—6.44%,  5.16%)  does  so  for  the  temperature  difference  between 
ends  of  the  modules. 

Furthermore,  the  computational  model  simulates  with  enough 
accuracy  the  performance  of  TSC  systems  under  transient  condi¬ 
tions.  This  fact  indicates  that  the  model  is  a  powerful  tool  to  eval¬ 
uate  the  response  of  the  system  under  changing  working 
conditions,  including  those  when  the  modules  do  not  generate 
enough  electric  power  to  operate  the  fan. 

In  conclusion,  the  performance  of  TSC  systems  was  assessed  in  a 
previous  paper;  now,  the  present  paper  goes  further  in  this 
research  line  and  presents  a  computational  tool  for  design  and 
development.  Therefore,  real  applications  of  this  technology  are 
expected  in  the  short  term.  Applications  intended  to  replace 
widespread  forced-convection  systems,  saving  significant  amounts 
of  electricity  and  reducing  the  corresponding  greenhouse  emis¬ 
sions,  in  line  with  the  global  fight  against  climate  change. 
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Appendix  A 

This  appendix  presents  the  procedure  devised  to  obtain  the 
thermal  resistance  of  the  dissipator  as  a  function  of  the  voltage 
supplied  by  the  thermoelectric  modules  to  the  fan.  As  explained  in 
Section  2.3.4,  this  thermal  resistance  includes  the  thermal  contact 
between  the  dissipator  and  the  cold  extender,  the  heat  conduction 
through  the  dissipator  and  the  heat  convection  between  the  dis¬ 
sipator  and  the  environment.  This  last  term  depends  directly  on  the 
air  flow  forced  through  the  fins  of  the  dissipator,  which  in  turn 
depends  on  the  electric  power  supplied  by  the  modules  to  the  fan. 
Therefore,  to  obtain  such  an  expression,  one  must  deal  with  both 
heat  transfer  and  fluid  mechanics  phenomena,  with  no  accurate 
analytical  solution  [27,28],  which  require  the  use  of  more  complex 
procedures,  such  as  the  one  presented  here. 

A.1.  Methodology 

Once  defined  the  characteristics  of  the  cold  extender  (di¬ 
mensions  and  materials),  the  characteristics  of  the  dissipator,  the 
contact  between  them,  and  the  type  of  fan  and  its  location  (all  of 
them  indicated  in  Sections  3.1.3  and  3.1.4),  one  finds  that  the 
thermal  resistance  of  the  dissipator  depends  only  on  the  volume 
flow  rate  of  the  air  forced  through  the  fins  [27],  Fluid-dynamics 
software  provides  the  cited  relation,  as  explained  in  A.2. 

This  volume  flow  rate  can  be  easily  deduced,  taking  into  account 
that  the  pressure  drop  of  the  air  flowing  through  the  fins  equals  the 
increase  in  pressure  provided  by  the  fan  [28], 

On  one  hand,  this  pressure  drop  depends  on  the  volume  flow 
rate,  and  the  relation  can  be  calculated  by  fluid-dynamics  software, 
as  A.3  details. 

On  the  other  hand,  the  increase  in  the  air  pressure  caused  by  the 
fan  is  also  a  function  of  the  volume  flow  rate,  and  their  relation  is 
provided  by  the  fan  manufacturer.  However,  this  information  is 
usually  available  only  for  the  case  when  the  fan  is  supplied  with 
maximum  voltage.  This  is  the  reason  why  A.4  shows  a  procedure  to 
obtain  the  cited  relation  for  different  voltages. 


Fig.  At..  Temperature  distribution  in  the  dissipator,  provided  by  Fluent-Ansys,  for 
37.39  m3/h  of  air  volume  flow  rate  (3  m/s). 


Table  A.l.  Thermal  resistance  of  the  dissipator  and  pressure  drop  for 
different  air  volume  flow  rates. 


vWJh)  Tcs-WK)  RcS(K/W)  A“Pr  (Pa) 


0  0  57.0 

0.5  6.23  20.9 

1  12.46  13.3 

2  24.93  9.5 

3  37.39  8.1 


1.20  0 
0.44  1.18 

0.28  3.74 

0.20  13.10 

0.17  75.91 


A.2.  Thermal  resistance  of  the  dissipator  as  a  function  of  air  volume 
flow  rate 

As  can  be  seen  in  the  TSC  prototype  presented  in  Fig.  1,  the  air 
stream  enters  the  fan  from  above,  flows  through  the  fins  and  the 
wind  tunnel,  and  leaves  the  dissipator  laterally.  To  simulate  these 
phenomena,  Fluent-Ansys  is  the  fluid-dynamics  software  selected. 
To  this  end,  the  finned  dissipator  and  the  air  volume  between  the 
dissipator  and  the  wind  tunnel  are  drawn,  meshed  and  simulated 
with  the  following  boundary  conditions: 

•  Environment  temperature  equal  to  298  K  (25  °C) 

•  The  air  stream  leaves  the  fan  at  environment  temperature  and 
fixed  velocity.  Air  velocities  of  0.5, 1, 2  and  3  m/s  are  simulated. 
Then,  the  volume  flow  rate  is  easily  derived,  given  that  the  fan 
cross-area  is  3.462  x  1(T3  m2. 

•  Mass  conservation  at  the  outflow  area  is  imposed. 

•  Over  the  base  surface,  where  the  dissipator  contacts  the  cold 
extender,  47.5  W  of  heat  transfer  rate  is  imposed,  representing 
the  heat  transferred  from  the  extender  to  the  dissipator.  (Qcoid ) 

Once  a  simulation  is  complete,  Fluent-Ansys  provides  the  tem¬ 
perature  distribution  in  the  dissipator.  Likewise,  it  provides  the 
average  temperature  on  the  surface  in  contact  with  the  cold  extender, 
required  to  calculate  the  thermal  resistance  of  the  dissipator  with  Eq. 
(A.1 ).  As  an  example,  Fig.  A.1  shows  the  temperature  distribution  in 
the  dissipator  for  3  m/s  of  air  velocity,  where  the  hottest  area  corre¬ 
sponds  obviously  to  the  surface  in  contact  with  the  cold  extender. 
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Rcs  Tem-  (A.1) 

ficold 

Finally,  Table  A.1  presents  the  thermal  resistance  of  the  dis- 
sipator  for  different  values  of  the  volume  flow  rate,  and  the  cor¬ 
relation  between  them  is  provided  by  Eq.  (A.2).  This  table  also 
includes  the  case  when  the  fan  is  not  powered,  so  no  air  stream  is 
forced  through  the  fins,  whose  corresponding  thermal  resistance 
(1.2  K/W)  was  calculated  experimentally.  The  inverse  of  this 
expression  is  represented  by  the  broken  curve  in  Fig.  A.2. 

Rcs  =  l/(~  0.0030V2  +  0.2447V  +  0.8330)  (A.2) 


Ra-‘  (W/K) 


V  (mVh) 


Fig.  A.2..  Calculation  of  the  thermal  resistance  of  the  dissipator  for  any  voltage  sup- 


air  volume  flow  rate,  for  12  V  supplied  to  the 


A+Pr12(Pa)  44.79  44.05  29.11  18.66  3.73  0 

Vu  (m3/h)  0  6.80  27.18  47.57  67.96  71.36 


A3.  Pressure  drop  as  a  function  of  the  volume  flow  rate 

The  pressure  drop  of  the  air  flowing  through  the  fins  is  calcu¬ 
lated  by  Eq.  (A.3),  composed  of  three  terms,  from  left  to  right: 
pressure  drop  of  the  inflow,  pressure  drop  due  to  friction  with  the 
fins  and  pressure  drop  of  the  outflow  [28],  Pressure  coefficients  for 
inflow  and  outflow  are  respectively  0.5  and  1;  air  density  is 
considered  1.254  kg/m3;  and  Fluent-Ansys  analysis  provides  the 
pressure  drop  due  to  friction. 

A  Pr  =  Kinfdf  +  A  PrMction  +  Kmtf^  (A.3) 

The  pressure  drop  is  calculated  for  different  air  velocities  and 
the  results  are  presented  in  Table  A.1.  In  this  table,  the  first  point  has 
been  added,  indicating  obviously  that  no  pressure  drop  occurs 
when  the  fan  is  not  powered.  Finally,  Eq.  (A.4)  shows  the  correla¬ 
tion  between  the  air  pressure  drop  and  its  volume  flow  rate,  which 
is  represented  by  the  solid  curve  in  Fig.  A.2. 

A  Pr  =  0.0001V4  -  0.0057V3  +  0.0890V2  +  0.1748V  (A.4) 

Table  A.3.  Experimental  values  of  angular  velocity  of  the  fan,  depending  on  the  sup¬ 
plied  voltage. 


V(V)  0  4  6  8  10  12 

to  (rpm)  0  950  1530  2040  2490  2850 


A.4.  Air  pressure  increase  versus  air  volume  flow  rate,  for  different 
voltages  supplied  to  the  fan 


The  fan  manufacturer  provides  a  curve  of  the  increase  in  the  air 
pressure  versus  its  volume  flow  rate,  for  maximum  voltage  sup¬ 
plied  to  the  fan  (12  V  for  the  fan  used  in  the  TSC  prototype).  To 
represent  this  curve,  we  have  extracted  the  information  presented 
in  Table  A.2,  which  leads  to  the  fitting  expression  given  by  Eq.  (A.5) 
and  represented  by  the  upper  dotted  curve  in  Fig.  A.2. 

A+Pr12  =  1.5  x  10-7Vj2  -  3.3  x  lfT5V?2  +  0.260647V^2 
-  0.08295V^2  +  0.3437V12  +  44.79 

(A.5) 

To  calculate  this  correlation  for  any  other  voltage,  one  must  use 
similarity  techniques  for  incompressible  fluids  [28],  These  indicate 
that  the  correlation  between  air  pressure  increase  provided  by  a  fan 
supplied  with  12  V,  and  the  air  pressure  increase  provided  by  a  fan 
with  any  other  voltage  can  be  determined  by  Eq.  (A.6).  Likewise,  Eq. 
(A.7)  provides  a  similar  expression  for  the  air  volume  flow  rate.  In 
both  expressions,  the  key  parameter  is  the  correlation  between 
angular  velocities  of  the  fan  at  different  voltages. 


A+Pr 

A+Pt12 


V 

Vl2 


Wl2 


(A.6) 

(A.7) 


Table  A.3  provides  the  angular  velocity  at  different  voltages,  all  of 
them  obtained  experimentally,  which  lead  to  the  fitting  expression 
presented  by  Eq.  (A.8).  It  is  important  to  note  that  the  fan  requires  at 
least  4  V  to  rotate  steadily,  so  Eq.  (A.8)  applies  only  to  these  cases.  For 
the  rest  ( V  <  4),  the  fan  is  considered  to  exhibit  null  angular  velocity. 


w  =  -8.9286V2  +  380.86V  -  432 


(A.8) 


■■etaLli 


V(V)  0  4  6  8  10  12 

«cs(  K/W)  1,2  0.45  0.22  0.20  0.19  0.18 


A.5.  Final  step:  thermal  resistance  of  the  dissipator  as  a  function  of 
the  voltage  supplied  to  the  fan 

Once  the  previous  correlations  have  been  obtained,  one  can 
derive  the  thermal  resistance  of  the  dissipator  for  a  specific  voltage 
with  a  procedure  similar  to  that  presented  in  Fig.  A.2,  which  dis¬ 
plays  the  cases  when  the  fan  is  powered  with  12  V  and  6  V.  The 
procedure  remains  as  follows: 

Firstly,  given  that  the  fan  is  supplied  with  6  V,  one  derives  the 
corresponding  angular  velocity  with  Eq.  (A.8).  Then,  Eqs.  (A.6)  and 
(A.7)  allow  us  to  obtain  a  new  table  similar  to  Table  A.2,  which  will 
display  the  increase  in  the  air  pressure  as  a  function  of  the  volume 
flow  rate,  for  6  V  of  supplied  voltage.  This  information  leads  to  the 
corresponding  fitting  expression,  represented  by  the  lower  dotted 
curve  in  Fig.  A.2.  Then,  the  point  where  this  new  curve  intersects 
the  curve  of  the  pressure  drop  provides  the  actual  volume  flow  rate. 
Finally,  Eq.  (A.2)  provides  the  thermal  resistance  corresponding  to 
this  volume  flow  rate,  as  the  broken  line  in  Fig.  A.2  displays. 

This  procedure  must  be  repeated  for  several  voltages  to  obtain 
Table  A.4,  from  which  one  can  easily  derive  the  final  expression  of 
the  thermal  resistance  of  the  dissipator  as  a  function  of  the  voltage 
supplied  to  the  fan,  finally  given  by  Eq.  (A.9). 

If  V>4V 

Rcs  =  V(-  0.0038V4  +  0.1412 V3  - 1 ,9399V2  + 1 1 .805 V-  22.023)  (K/W) 

If  V<4V 
Rcs  =  1.2(K/W) 

(A.9) 
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